Based on the replication code for this paper, this workbook aims to explore Reviewer 2s objection:

“All legislators that remain in the dataset between t and t+1 increase in experience 1 unit. Thus, there is no identifying variation to separate the effect of experience from the effect of time shocks in the within-individual specification.

Thus, the coefficients on the dummy variables indicating different levels of legislative experience are not identified in the OLS regression with time- and legislator fixed effects.”

Replication data

From code/replication.R, which starts with the main data file (dcounts_min.Rdata) and makes transformations for three sets of models:

counts per member-year (dcounts_tenure.Rdata/AgencyComm.dta)
ratio per member-year (dcounts_ratio.Rdata/ProportionContact.dta)
counts per district-year (dcounts_per_district.Rdata, DistrictLevel.dta)

The same three datasets are used in Replication.do

load(here::here("data", "dcounts_tenure.Rdata"))
load(here::here("data", "dcounts_ratio.Rdata"))
load(here::here("data", "dcounts_per_district.Rdata"))
load(here::here("data" , "members.Rdata"))


# var names for tables
dcounts_tenure %<>% 
  mutate(Legislator = icpsr,
         Legislator_x_Agency = icpsr_agency,
         Year_x_Agency = agency_year)

dcounts_tenure %>% head() %>% kablebox()

agency	icpsr	chamber.x	year	congress	party_code	state_abbrev	district_code	nominate.dim1	presidents_party	female	chair	ranking_minority	majority	prestige	prestige_chair	yearelected	state	first_cong	first_year	tenure	chamber.y	oversight_committee	oversight_committee_chair	icpsr_agency	agency_year	max_year	survive	chamber	n	Legislator	Legislator_x_Agency	Year_x_Agency
ABMC	10713	House	2016	114	100	MI	13	-0.658	1	0	0	1	0	0	0	2014	michigan	89	1965	51	NA	NA	NA	ABMC_10713	ABMC_2016	53	1	House	1	10713	ABMC_10713	ABMC_2016
ABMC	10713	House	2017	115	100	MI	13	-0.658	0	0	0	1	0	0	0	2016	michigan	89	1965	52	NA	NA	NA	ABMC_10713	ABMC_2017	53	1	House	1	10713	ABMC_10713	ABMC_2017
ABMC	10713	House	2018	115	100	MI	13	-0.658	0	0	0	1	0	0	0	2016	michigan	89	1965	53	NA	NA	NA	ABMC_10713	ABMC_2018	53	1	House	1	10713	ABMC_10713	ABMC_2018
ABMC	13035	House	2016	114	100	NY	13	-0.514	1	0	0	0	0	1	0	2014	new york	92	1971	45	NA	NA	NA	ABMC_13035	ABMC_2016	45	1	House	1	13035	ABMC_13035	ABMC_2016
ABMC	14009	Senate	2016	114	200	MS	0	0.287	0	NA	1	0	1	1	1	2014	mississippi	93	1973	43	NA	NA	NA	ABMC_14009	ABMC_2016	45	1	Senate	1	14009	ABMC_14009	ABMC_2016
ABMC	14009	Senate	2017	115	200	MS	0	0.287	1	NA	1	0	1	1	1	2016	mississippi	93	1973	44	NA	NA	NA	ABMC_14009	ABMC_2017	45	1	Senate	1	14009	ABMC_14009	ABMC_2017

Reviewer 2s objection is about the models using (1) counts per member-year, so that is the DV I focus on here.

Values to predict at

The most frequent legislator agency pair: McCain x VA with 2159 letters.

The average is 73 letters per legislator per year, a little over 1.29 letters per agency per year.

For predictions, we need an average member-year-agency triad. The first table below shows the members who averaged 73, entered Congress between 2007 and 2011, were observed in our data for at least seven years, and most often submitted exactly one letter to an agency (the closest integer 1.29). For realism, this person would ideally, at some point in our data, become a chair. Of the modal legislators, Phil Roe meets this last criterion. The IRS is a fairly modal agency and one of the agencies to which Phil Roe submitted one letter in 2015, so predictions are at Phil Roe X IRS X 2015. Other legislators, years, and agencies are simply intercept shifts.

modal <- dcounts_tenure %>% 
  group_by(icpsr) %>% 
  filter(first_year > 2006, first_year < 2013, max_year > 6) %>% 
  mutate(Chair = sum(chair, na.rm = T) >1) %>% 
  group_by(icpsr, year) %>%
  mutate(perLeg = sum(perYear)) %>% 
  filter(perYear == 1, perLeg == 73 ) %>% 
  select(perAgency = perYear, perLeg, icpsr, icpsr_agency, agency_year, first_year, max_year, Chair) %>% 
  distinct() %>% 
  ungroup() %>% left_join(members %>% select(bioname, icpsr))



modal %<>% filter(Chair)


# by chair and presidents' party
values <- tidyr::expand(dcounts_tenure,
                        agency = "Treasury_IRS",
                        chair = chair, 
                        ranking_minority    = FALSE, 
                        prestige = TRUE, 
                        first = FALSE,
                        second = FALSE,
                        third = FALSE,
                        fourth = FALSE,
                        fifth = FALSE, 
                        sixth = FALSE,
                        Legislator_x_Agency = "Treasury_IRS_20947",
                        Year_x_Agency = "Treasury_IRS_2015",
                        Legislator = "20947",
                        Year = "2015",
                        majority = TRUE,
                        presidents_party = presidents_party
                   ) %>% 
  drop_na(majority, chair)

# by year 
values_tenure <- tidyr::expand(dcounts_tenure,
                        agency = "Treasury_IRS",
                        chair = chair, 
                        ranking_minority    = FALSE, 
                        prestige = TRUE, 
                        first = first,
                        second = second,
                        third = third,
                        fourth = fourth,
                        fifth = fifth, 
                        sixth = sixth,
                        Legislator_x_Agency = "Treasury_IRS_20947",
                        Year_x_Agency = "Treasury_IRS_2015",
                        Legislator = "20947",
                        Year = "2015",
                        majority = TRUE,
                        presidents_party = TRUE
                   ) %>% 
  drop_na(chair)

Simulated data

Counts per legislator per agency per year

I replace the dependent variable, the count of letters that each legislator writes to each agency each year with a normal distribution with a mean 1 (the true mean in the data is 1.29).

I also introduce a correlation with the calendar year by adding +1 for each year beyond 2007 (the first year in the data.) In this simulation 2017 would have a mean of 11.

Later, I will introduce a correlation with tenure. (See “100 simulations section” below)

n <- nrow(dcounts_tenure)

sim <- dcounts_tenure %>% 
  # REPLACE DV WITH NORMAL DIST
  mutate(perYear = rnorm(n, year-2007)
         # OTHER VARS LEFT ALONE 
         # chair,
         # first,
         # second,
         # third,
         # fourth,
         # fifth,
         # sixth,
         # Year, 
         # Agency,
         # Legislator
  ) 

dcounts_tenure <- sim

Models

Our primary models are a series of difference-in-differences regressions that examine changes that are within legislator and agency pairs.

Estimates from Model 2 (Column 2 below) provides the estimated effects from the difference-in-differences specification.

$Y_{ijt} = {}^{’} {it} + {s = 1}^{6} {s} ({it} = s) + {ij} + {jt} + m_{it} + p_{it} + _{ijt} $

Where $Y_{ijt}$ represents the number of requests legislator $i$ makes to agency $j$ in year $t$. Our analysis in this section is at the legislator-agency-year level. $\gamma_{ij}$ is a fixed effect for the legislator-agency pair. This fixed effect accounts for legislators’ characteristics, such as legislators who are more skillful at filling constituency service requests than other legislators. Critically for our research design, this fixed effect also enables us to account for time-invariant constituent demands, ensuring differences in constituent demand do not drive our results. It also accounts for state and districts characteristics, including population, demographics, and local industries that might be particularly likely to request help with specific agencies. This difference-in-difference design ensures that coefficients ${\beta}$ capture variation related to changes in institutional power or experience, not other factors that may vary across districts, legislators, or agencies. The model also accounts for the different periods for which data were available from each agency. $\delta_{jt}$ is an agency-year fixed effect. This takes into account agency-level shocks that may affect legislator requests.

Assuming that legislators’ trends in the level of requests follow parallel paths, $\boldsymbol{\beta}$ represents the average effect of changing institutional power on a legislator’s provision of constituency service.

In this same regression we also include indicators for legislators’ first six years in Congress, $ {s = 1}^{6} {s} _{it}$. The effects of interest $\eta_{1}, \eta_{2} : \eta_{6}$ describe how a legislator’s provision of constituency service at levels of seniority between one and six years differ from legislators who serve beyond six years. We focus on constituency service levels in each of the first six years of a legislator’s tenure to assess how constituency service changes over their initial years in Congress. This design allows us to assess the extent to which new legislators face start-up costs.

Estimation

Fixed effects general linear model fixest::fixest

feglm (perYear ~ 
 chair + ranking_minority + prestige +
 first + second + third + fourth + fifth + sixth + 
 majority + presidents_party | Year_x_Agency + Legislator_x_Agency, 
 cluster = "Legislator",
 data = dcounts_tenure)

Where

perYear = $Y_{ijt}$ = Letters per legislator $i$ per agency $j$ per year $t$

chair= $\textbf{Committee Position}_{it}$ = Legislator $i$ is a committee chair in year $t$

first + second + third + fourth + fifth + sixth = $\text{I}\left(\text{tenure}_{it} = s\right)$ = Legislator i is in their $s$ (first, second, third, fourth, fifth, sixth) year of service in calendar year $t$

$\gamma_{ij}$ = a fixed effect for the legislator-agency pair = Legislator_x_Agency

$\delta_{jt}$ = an agency-year fixed effect = Year_x_Agency

Counts per member per year

Clustering standard errors at the legislator level:

“Clustering on the panel variable produces an estimator of the VCE that is robust to cross-sectional heteroskedasticity and within-panel (serial) correlation that is asymptotically equivalent to that proposed by Arellano (1987)”

Replicated using fixest::feglm here and xtreg in Replication.do. See more about robust (clustered) standard errors with fixed effects are calculated by fixest here and by xtreg here.

# Coef Map
cm = c("chair" = "Committee Chair",
       "ranking_minority" = "Ranking Member",
       "prestige" = "Prestige Committee",
       "first" = "First Year",
       "second" = "Second Year",
       "third" = "Third Year",
       "fourth" = "Fourth Year",
       "fifth" = "Fifth Year",
       "sixth" = "Sixth Year",
       "majority" = "Majority",
       "presidents_party" = "President's Party",
       "Legislator" = "Legislator", 
       "Agency" = "Agency",
       "FE: icpsr_agency" = "Legislator × Agency Fixed Effects",
       "icpsr_agency" = "Legislator × Agency Fixed Effects")

coef_omit = "(Intercept)|majority|presidents_party"
coef_omit = "none"

setFixest_dict(cm)

gof_omit = "R.*|AIC|BIC|Log.*|Std.*"

# table formatting to match stata (sort of)
format_table <- . %>% 
  row_spec(row = 1, bold = T, hline_after = TRUE) %>% 
      kable_styling(font_size = 11) %>% 
                    #full_width = TRUE, 
                    #latex_options = c("repeat_header")) %>%
  str_replace("Num.Obs.", "Observations") %>% 
  str_replace("Std.Errors", "Robust/Clustered Std. Errors") %>% 
  str_replace("FE: Legislator_x_Agency", "Legislator x Agency FE") %>%
  str_replace("FE: Year_x_Agency", "Year x Agency FE") %>% 
  str_replace("FE: Year", "Year Fixed Effects") %>% 
  str_replace("FE: Legislator", "Legislator Fixed Effects") %>%
    str_replace_all("X|✓", "\\\\\\checkmark") %>% 
  # a random midrule appeared in the wrong place 
  #str_remove("\\midrule") %>% 
  #  extract just the table, no caption etc
  str_extract("\\\\begin\\{tabular\\}(\n|.)*tabular\\}")

Total Letters

Coeficient Plots

Figures: figs/m-total-[1:4].png

“cross-sectional”
“diff-in-diff”
“cross-sectional re-elected”
“diff-in-diff re-elected”

testing = F

# paper table 
# 1 
# cross-sectional 
m_total_cross <-feglm (perYear ~ 
                         chair + ranking_minority + prestige + 
                    first + second + third + fourth + fifth + sixth + 
                    majority + presidents_party | Year_x_Agency,
                    cluster = "Legislator", 
           data = dcounts_tenure)

if(testing){modelsummary(m_total_cross, gof_omit = "R.*", coef_omit = coef_omit)}

coefplot(m_total_cross, horiz = T, drop = "(Intercept)") 



# 2
m_total_dnd <-feglm (perYear ~ 
                       chair + ranking_minority + prestige + 
                    first + second + third + fourth + fifth + sixth + 
                    majority + presidents_party | Year_x_Agency + Legislator_x_Agency, 
                    cluster = "Legislator",
                    # ALT vcov = hetero ~ ssc(cluster.adj = TRUE),
           data = dcounts_tenure)

if(testing){modelsummary(m_total_dnd, gof_omit = "R.*", coef_omit = coef_omit)}

coefplot(m_total_dnd, horiz = T) 


# 3
m_total_2nd <-feglm (perYear ~ 
                       chair + ranking_minority + prestige + 
                    first + second + third + fourth + fifth + sixth + 
                    majority + presidents_party | Year_x_Agency + Legislator_x_Agency, 
                    cluster = "Legislator", 
           data = dcounts_tenure %>% filter(survive == 1))

if(testing){modelsummary(m_total_2nd, gof_omit = "R.*", coef_omit = coef_omit)}
coefplot(m_total_2nd, horiz = T) 


# 4
m_logtotal_dnd <-feglm (log(perYear + 1) ~ 
                          chair + ranking_minority + prestige + 
                    first + second + third + fourth + fifth + sixth + 
                    majority + presidents_party | Year_x_Agency + Legislator_x_Agency, 
                    cluster = "Legislator", 
           data = dcounts_tenure)

if(testing){modelsummary(m_logtotal_dnd, gof_omit = "R.*", coef_omit = coef_omit)}

coefplot(m_logtotal_dnd, horiz = T)

Total Letters Per Year

Total Letters Table

Table: tables/models_total.tex (sans stars)

models_total <- list(
  "(1)" = m_total_cross,
  "(2)" = m_total_dnd,
  "(3)" =  m_total_2nd,
  "(4)" = m_logtotal_dnd
)

rows <- tibble(
  term = c("Dependent Variable", 
           "Majority",
           "President's Party",
           "All Legislators", 
           "Served At Least 2nd Term"),
  `(1)` = c("Count", "✓", "✓", "✓", ""),
  `(2)` =c("Count", "✓", "✓", "✓", ""),
  `(3)` = c("Count", "✓","✓","", "✓"),
  `(4)` = c("Log(Count+1)","✓","✓", "✓", "") 
)

rows <- tibble(
  term = c("Dependent Variable", 
           #"Majority",
           #"President's Party",
           "All Legislators", 
           "Served At Least 2nd Term"),
  `(1)` = c("Count", "✓", ""),
  `(2)` =c("Count", "✓", ""),
  `(3)` = c("Count","", "✓"),
  `(4)` = c("Log(Count+1)", "✓", "") 
)

# subset
attr(rows, 'position') <- c(0, 20,21,22,23)

# call coefs
attr(rows, 'position') <- c(0, 24,25)


# HTML
modelsummary(models_total, stars = T, 
             add_rows = rows, 
             coef_map = cm,
             gof_omit = gof_omit,
             coef_omit = coef_omit,
             title = "The Effect of Expierence and Institutional Power on Legislator Contacts", 
             notes = list("Robust standard errors in parentheses, clustered by legislator.")) %>% row_spec(row = 1, bold = T, hline_after = TRUE)

The Effect of Expierence and Institutional Power on Legislator Contacts
	(1)	(2)	(3)	(4)
Dependent Variable	Count	Count	Count	Log(Count+1)
Committee Chair	0.003	0.014	0.013	0.004
	(0.006)	(0.010)	(0.010)	(0.003)
Ranking Member	−0.008	−0.006	−0.006	−0.002
	(0.006)	(0.009)	(0.009)	(0.003)
Prestige Committee	−0.001	0.009	0.008	0.000
	(0.003)	(0.007)	(0.007)	(0.002)
First Year	−0.001	−0.005	−0.004	0.000
	(0.006)	(0.011)	(0.011)	(0.005)
Second Year	−0.006	−0.009	−0.008	−0.003
	(0.006)	(0.011)	(0.011)	(0.004)
Third Year	0.007	0.007	0.007	0.004
	(0.007)	(0.010)	(0.010)	(0.004)
Fourth Year	−0.004	−0.004	−0.005	−0.002
	(0.007)	(0.010)	(0.010)	(0.003)
Fifth Year	0.014+	0.014	0.015	0.003
	(0.007)	(0.009)	(0.010)	(0.003)
Sixth Year	0.000	−0.001	−0.001	−0.001
	(0.007)	(0.009)	(0.009)	(0.003)
Majority	−0.003	−0.005	−0.005	−0.004*
	(0.003)	(0.005)	(0.005)	(0.002)
President’s Party	−0.001	−0.001	−0.002	0.002
	(0.003)	(0.004)	(0.004)	(0.002)
All Legislators	✓	✓		✓
Served At Least 2nd Term			✓
Num.Obs.	412111	412111	388997	407272
FE: Year_x_Agency	X	X	X	X
FE: Legislator_x_Agency		X	X	X
+ p < 0.1, * p < 0.05, p < 0.01, * p < 0.001
Robust standard errors in parentheses, clustered by legislator.

# TeX 
modelsummary(models_total, # stars = T, 
             add_rows = rows,               
             coef_map = cm, 
             gof_omit = gof_omit, 
             coef_omit = coef_omit,
             output = "latex",
             #title = "\\\\footnotesize The Effect of Expierence and Institutional Power on Legislator Contacts", 
             notes = list("\\\\footnotesize Robust standard errors in parentheses, clustered by legislator.")) %>% 
  format_table() %>%
  write_lines(file = here::here("docs", "tables", "models_total.tex") )


beta <- m_total_dnd$coefficients %>%
  round(3) %>% 
  as_tibble(rownames = "beta") %>% 
  pivot_wider(names_from = beta)

se <- m_total_dnd$se %>%
  round(3) %>% as_tibble(rownames = "se") %>%
  pivot_wider(names_from = se)

Total Letters Predictions

Figures: figs/m-total-predicted-[1:4].png

# helper functions to plot predicted values
ggchair <- function(predicted = predicted) {
  
  predicted$chair %<>% str_replace("0", "Not Chair") %>% str_replace("1", "Chair")
  
  predicted$presidents_party %<>% str_replace("0", "Not President's Party") %>% str_replace("1", "President's Party")

predicted %>%  
  ggplot() + 
  aes(x = "", 
      y = estimate, 
      fill = chair,
      shape = chair,
      color = chair,
      ymin = conf.low,
      ymax = conf.high) + 
  geom_pointrange(position =  position_dodge(width = -.5)  )  + 
  scale_fill_viridis_d(begin = 0, end = .6, option = "cividis") +
  scale_color_viridis_d(begin = 0, end = .6, option = "cividis") +
      theme_minimal() + 
  theme(panel.border  = element_blank(),
        panel.grid.major.x = element_blank()) + facet_grid(. ~ presidents_party) 
}


ggtenure <- function(predicted = predicted) {
  
predicted %<>%
  # drop estimates from impossible values
      group_by(rowid, estimate, conf.low, conf.high, chair) %>% 
  mutate(sum = sum(first, second, third, fourth, fifth, sixth),
         more = ifelse(sum == 0, 1,0)) %>% 
    filter(sum < 2) %>% 
    pivot_longer(cols = c("first", "second", "third", "fourth", "fifth", "sixth", "more")) %>% 
    select(name, value) %>% 
    filter(value == 1) %>% 
    # clean up for presentation
    mutate(year_in_congress = name %>% 
             str_replace("more", "7\nor more") %>% 
             str_replace("sixth", "6") %>% 
             str_replace("fifth", "5") %>% 
             str_replace("fourth", "4") %>% 
             str_replace("third", "3") %>% 
             str_replace("second", "2") %>% 
             str_replace("first", "1") 
             ) %>% 
    ungroup()
  
# Ideally, we could plot against data, but there is so much variation that you can no longer distinguish differences in predicted values 
# predicted %<>% full_join(dcounts_tenure2 %>% mutate(predicted = NA))
  
predicted %<>% filter(chair == 0 | name %in% c("sixth", "more"))

  predicted$chair %<>% str_replace("0", "Not Chair") %>% str_replace("1", "Chair")
  

predicted %>%  
  ungroup() %>% 
  ggplot() + 
  aes(x = year_in_congress, 
      y = estimate, 
      shape = chair,
      color = chair,
      ymin = conf.low,
      ymax = conf.high) + 
  geom_pointrange(position =  position_dodge(width = -.3)  )  + 
  scale_fill_viridis_d(begin = 0, end = .6, option = "cividis") +
  scale_color_viridis_d(begin = 0, end = .6, option = "cividis") +
      theme_minimal() + 
  theme(panel.border  = element_blank(),
        panel.grid.major.x = element_blank()) 
  
}

predicted <- predictions(m_total_cross,
                         newdata = values)

# Cross-sectional predictions
predicted %>%
    mutate(estimate = estimate*90,
         conf.low = conf.low*90, 
           conf.high = conf.high*90) %>%
  ggchair() + 
  labs(title = "Predicted Total Letters per Year\n(Cross-Sectional)",
       x = "",
       y = "Predicted Total Letters Per Year",
       fill = "",
       color = "",
       shape = "") 

# by tenure
predicted <- predictions(m_total_cross,
                         newdata = values_tenure)

predicted %>%
    mutate(estimate = estimate*90,
         conf.low = conf.low*90, 
           conf.high = conf.high*90) %>% 
ggtenure() +
  labs(title = "Predicted Total Letters Per Year\n(Cross-Sectional)",
       x = "Years Serving in Congress",
       y = "Predicted Letters Per Year",
       fill = "",
       color = "",
       shape = "") 


# diff in diff 
predicted <- predictions(m_total_dnd,
                         newdata = values)


# predictions by chair and presidents' party
predicted %>%
    mutate(estimate = estimate*90,
           conf.low = conf.low*90, 
           conf.high = conf.high*90) %>%
  ggchair() + 
  labs(title = "Predicted Total Letters per Year\nDifference in Differences (Within-Legislator)",
       x = "",
       y = "Predicted Letters Per Year",
       fill = "",
       color = "",
       shape = "") 

# by tenure
predicted <- predictions(m_total_dnd,
                         newdata = values_tenure)

predicted %>%
    mutate(estimate = estimate*90,
           conf.low = conf.low*90, 
           conf.high = conf.high*90) %>%
ggtenure() +
    geom_line(aes(x = as.numeric(year_in_congress %>% str_sub(1,1)))) + 
  labs(title = "Predicted Total Letters Per Year\nDifference in Differences (Within-Legislator)",
       x = "Years Serving in Congress",
       y = "Predicted Letters Per Year",
       fill = "",
       color = "",
       shape = "")

Total

100 simulations

Subletting data to legislators in their first 6 years and setting the mean of the simulated DV (letters per legislator per year per agency) to their year (tenure) in office.

# A function to estimate diff in diff with a new random sample for DV 
sim_data <- dcounts_tenure %>% filter(tenure<7, tenure>0)

# sim_data$tenure

Cross sectional

simulation_cross <- function(x){
  
  set.seed(x)
  
  data <- sim_data %>% 
  # REPLACE DV WITH NORMAL DIST
  mutate(perYear = rnorm(n, tenure)) # year-2007))

m_total <-feglm (perYear ~ 
                       #chair + ranking_minority + prestige + 
                    first + second + third + fourth + fifth + sixth | Year_x_Agency, 
                    cluster = "Legislator",
                    # ALT vcov = hetero ~ ssc(cluster.adj = TRUE),
           data = data)

beta <- m_total$coefficients %>%
  #round(3) %>% 
  as_tibble(rownames = "beta") %>% 
  pivot_wider(names_from = beta)

return(beta)
}

# run 100 times 
sims <- map_dfr(.x = 1:100, .f =  simulation_cross)

# # plot estimates
# sims %>% 
#   ggplot() +
#   aes(x = first) + 
#   geom_density(fill = "grey")  #+  geom_vline(xintercept = 1)

sims %>% 
  ggplot() +
  aes(x = second) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = third) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = fourth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = fifth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = sixth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

Diff in diff

simulation_dnd <- function(x){
  
  set.seed(x)
  
  data <- sim_data %>% 
  # REPLACE DV WITH NORMAL DIST
  mutate(perYear = rnorm(n, tenure)) # year-2007))

m_total_dnd <-feglm (perYear ~ 
                      # chair + ranking_minority + prestige + 
                      # majority + presidents_party 
                      first + second + third + fourth + fifth + sixth | Year_x_Agency + Legislator_x_Agency, 
                    cluster = "Legislator",
                    # ALT vcov = hetero ~ ssc(cluster.adj = TRUE),
           data = data)

beta <- m_total_dnd$coefficients %>%
  #round(3) %>% 
  as_tibble(rownames = "beta") %>% 
  pivot_wider(names_from = beta)

return(beta)
}

# run 100 times 
sims <- map_dfr(.x = 1:100, .f =  simulation_dnd)

# plot estimates
# sims %>% 
#   ggplot() +
#   aes(x = first) + 
#   geom_density(fill = "grey")  #+  geom_vline(xintercept = 1)

sims %>% 
  ggplot() +
  aes(x = second) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = third) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = fourth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = fifth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

sims %>% 
  ggplot() +
  aes(x = sixth) + 
  geom_density(fill = "grey") #+ geom_vline(xintercept = 2)

Reproduceable Example Code